%%%%%%%%%%%%
% Pricing regressions
% Pricing regression, and return factor model regression
% Saves various outputs for analysis (performance, market price, portfolio weights, return)
% Liuren Wu, liurenwu@gmail.com
%%%%%%%%%%%%
format compact; clearvars;
load('./data/Momentsw.mat');
load('./data/PLw.mat','PLD');

tic;
%BMS greeks, with F normalized to 100
[Theta, dVega, dGamma, dVanna, dVolga] = FunForwardVGVVGreeks(100, d1w, Iw, matw);


d1m=repmat(d1(:)',nm,1);d1m=(d1m(:)); nd1m=normcdf(d1m);
Nd1m=normcdf(d1m);
mm=repmat(mat(:),1,nk);mm=(mm(:));smm=sqrt(mm);lmm=log(mm);
pdm=repmat(pdelta(:)'*100,nm,1); pdm=pdm(:);

nmk=nm*nk;%cross section size
bdsd=2;bdsm=1;

np=4;
R2=NaN(Tw,nh,nc);
B2V=NaN(Tw,np,nc);
VC=NaN(Tw,np+1,nc);

RRP=NaN(Tw,np,nc);NRP=RRP;
ev=NaN(Tw,nmk,nc);evv=ev;

CSMoments=NaN(Tw,nmk,3,nh+1,nc);
FPL=NaN(Tw,np+1,2,nc);
BO=NaN(Tw,nc);
NOBV=NaN(Tw,nh+1,nc);
FactorModel=NaN(Tw,np+2,nc);
RVC=NaN(Tw,np+1,nc);ARVC=RVC;
wr=NaN(nmk,np+1,Tw);
wx=NaN(nmk,np,Tw);


nz=4;
MSEOS=NaN(Tw,16,nc);
MSEOZ=NaN(Tw,16,nc);

AVGM=NaN(Tw,3,nc);
tradeflag=NaN(Tw,nc);

rng('default');
for t=T0+1:Tw
    %randomly generated in sample and out-of-sample index (for appendix
    %performance comparison)
    rdt=rand(nmk,1);
    idxis=rdt<=prctile(rdt,80); %in-sample index
    idxos=idxis==0; %out-of-sample index

    for c=1:nc
        y=reshape(-Theta(t,:,:,c),nmk,1);
        
        vegat=reshape(dVega(t,:,:,c),nmk,1);
        gammat=reshape(dGamma(t,:,:,c),nmk,1);
        vannat=reshape(dVanna(t,:,:,c),nmk,1);
        volgat=reshape(dVolga(t,:,:,c),nmk,1);
        It=reshape(Iw(t,:,:,c),nmk,1);

        xx=NaN(nmk,np,nh);
        mom=NaN(nmk,3,nh);
        rhov=NaN(nh,1);
        for h=1:nh
            hvt=reshape(HVw(t,:,:,c,h),nmk,1);
            hgt=reshape(HGw(t,:,:,c,h),nmk,1);
            hot=reshape(HOw(t,:,:,c,h),nmk,1);
            
            %cross-sectional smoothing on moment conditons
            hvt= FunMomentCSSmoothing(hvt,Nd1m,lmm,bdsd,bdsm);
            hot= FunMomentCSSmoothing(hot,Nd1m,lmm,bdsd,bdsm);
            hgt= FunMomentCSSmoothing(hgt,Nd1m,lmm,bdsd,bdsm);

            rhov(h)=nanmean(hgt(:))./sqrt(nanmean(hvt(:)).*nanmean(hot(:)));
            hovt=sqrt(hot.*hvt);

            mom(:,:,h)=[hvt,hgt,hot];
            X=[vegat.*(hot),  0.5*gammat.*hvt, 0.5*volgat.*hot,vannat.*(hovt)];
            ii=isfinite(sum(X,2)+y);
            Ni=sum(ii);
            NOBV(t,h,c)=Ni;
            if Ni>15
                yi=y(ii); Xi=X(ii,:);
                B=Xi\yi;
                ei=yi-Xi*B;  R2t=1-mean(ei.^2)/var(yi);
                R2(t,h,c)=R2t;
                xx(ii,:,h)=Xi;
            else
                %disp(datestr(wdate(t)))
                %Ni
                %pause
            end
        end
        [R2t,indmr]=max(R2(t,:,c)); w=zeros(3,1); w(indmr)=1;
        mom(:,:,4)=mom(:,:,1)*w(1)+mom(:,:,2)*w(2)+mom(:,:,3)*w(3);
        CSMoments(t,:,:,:,c)=mom;
        X=xx(:,:,1)*w(1)+xx(:,:,2)*w(2)+xx(:,:,3)*w(3); %regressors

        Bnull=[0,1,1,rhov(1)]';
        BO(t,c)=rhov(1);
        ii=isfinite(sum(X,2)+y);
        Ni=sum(ii);NOBV(t,h+1,c)=Ni;
        if Ni>15
            tradeflag(t,c)=1;
            AVGM(t,:,c)=mean(mom(ii,:,4));

            iis=ii&idxis; ios=ii&idxos;

            %full-sample estimation
            yi=y(ii);
            Xi=X(ii,:);
            B=Xi\yi;
            e=y-X*B;ev=e./gammat./It;
            ei=e(ii); eiv=ev(ii); %full sample
            eios=e(ios);eivos=ev(ios); %full sample estimates on out of sample

            %%in-sample out of sample
            yii=y(iis);
            Xii=X(iis,:);
            Bi=Xii\yii;
            e=y-X*Bi;
            evs=e./gammat./It;
            eis=e(iis); eos=e(ios);
            evis=evs(iis); evos=evs(ios);



            MSEOS(t,:,c)=[mean(ei.^2), mean(eis.^2), mean(eos.^2),mean(eios.^2), ...
                mean(yi.^2), mean(y(iis).^2), mean(y(ios).^2), mean(y(ios).^2), ...
                mean(eiv.^2),mean(evis.^2),mean(evos.^2),mean(eivos.^2), ...
                mean(abs(eiv)),mean(abs(evis)),mean(abs(evos)),mean(abs(eivos))];

            %%in-sample out-of-sample stats


            R2t=1-mean(ei.^2)/mean(yi.^2);
            eiv=ei./gammat(ii)./It(ii);

            wx(ii,:,t,c)=Xi/(Xi'*Xi); %portfolio weight
            cXX=(Xi'*Xi);cy=yi'*yi;
            aat= [((cXX*B/cy).*B);  R2t]; %variance contribution
            VC(t,:,c)=aat;

            B2V(t,:,c)=B; 
            ev(t,ii,c)=ei;   %theta/return space error
            evv(t,ii,c)=eiv; %vol space error

            RRP(t,:,c)=B-Bnull; %market price of risk
            Hr=[Xi, ei];              portwr=Hr/(Hr'*Hr); %portfolio weight
            swrt=sum(abs(portwr));    portwn=portwr./swrt; %fixed notional portfolio weight
            NRP(t,:,c)=RRP(t,:,c)./swrt(1:np); %market price per unit notional

            wr(ii,:,t,c)=portwr; 

            %% Alternative pricing model (appendix) under common risk assumption
            Z=[vegat, 0.5*gammat, 0.5*volgat,vannat];
            yi=y(ii);
            Zi=Z(ii,:);
            C=Zi\yi;
            ez=y-Z*C; ezv=ez./gammat./It;
            eiz=ez(ii);eizv=ezv(ii);
            eizos=ez(ios);eizvos=ezv(ios); %full sample estimates on out of sample

            %in-sample & out of sample
            yii=y(iis);
            Zii=Z(iis,:);
            Ci=Zii\yii;
            e=y-Z*Ci;
            evs=e./gammat./It;
            eis=e(iis); eos=e(ios);
            evis=evs(iis); evos=evs(ios);

            MSEOZ(t,:,c)=[mean(eiz.^2),mean(eis.^2), mean(eos.^2),mean(eizos.^2), ...
                mean(yi.^2), mean(y(iis).^2), mean(y(ios).^2), mean(y(ios).^2), ...
                mean(eizv.^2),mean(evis.^2),mean(evos.^2),mean(eizvos.^2), ...
                mean(abs(eizv)),mean(abs(evis)),mean(abs(evos)),mean(abs(eizvos))];


            %% return factor model
            PLt=reshape(PLD(t,:,:,c),nmk,1); %delta hedged return over next week
            PLi=PLt(ii);
            ij=isfinite(PLi);
            if sum(ij)>15
                %regression
                Hi=[Xi,ei];
                Hj=Hi(ij,:);yj=-PLi(ij);
                Bj=Hj\yj;
                ej=yj-Hj*Bj;R2j=1-mean(ej.^2)/mean(yj.^2);

                FactorModel(t,:,c)=[Bj;R2j];

                cXX=(Hj'*Hj);cy=yj'*yj;
                aat= [((cXX*Bj/cy).*Bj)];
                RVC(t,:,c)=aat;

                for k=1:np+1 %risk-adjusted factor portfolio
                    w0=portwr(:,k);
                    wn=portwn(:,k);
                    ww=[w0,wn];

                    FPL(t,k,:,c)=-PLi(ij)'*ww(ij,:); %PL for next day, recorded on today
                end

            end

        end

    end
end
toc;
save('./temp/PricingResultspartial.mat', ...
    'wdate','Tw','T0','VC','FPL','RRP','NRP','CSMoments','BO','ev','evv','FactorModel', 'RVC','wr','wx', ...
     'MSEOS', 'MSEOZ','AVGM','tradeflag',  ...
    'np','nmk','-mat');

return


